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Structural  Dynamics  Branch,  Structures  Division,  Flight  Dynamics  Directorate,  Wright 
Laboratory,  Wright-Patterson  Air  Force  Base,  Ohio.  The  study  reported  herein  was 
conducted  under  Work  Unit  2401TI00,  “Structural  Technology  Integration”. 

This  manuscript  was  released  in  September  1994  for  publication  as  a  Technical 
Memorandum  covering  work  performed  during  June  and  July  1994. 
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Summary 


This  Technical  Memorandum  presents  an  algorithm  for  deforming  a  two-  or  three- 
dimensional  aerodynamic  grid  given  deflections  on  a  physical  boundary.  This  algorithm  is 
to  be  used  in  coupling  a  computational  fluid  dynamics  code  with  a  computational 
structural  mechanics  code.  This  algorithm  is  intended  to  assist  engineers  in  the  Flight 
Dynamics  Directorate  with  the  analysis  of  aeroelastic  problems  using  state-of-the-art 
computational  methods.  The  deformation  algorithm  was  tested  successfully  on  both  two- 
and  three-dimensional  grids.  The  method  sufficiently  maintains  grid  quality  for  smooth 
deflections  of  realistic  orders  of  magnitude.  The  method  is  computationally  efficient:  the 
time  required  to  deform  the  grid  is  small  compared  to  the  time  required  to  solve  the  fluid 
dynamics. 
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Introduction 


The  Flight  Dynamics  Directorate  is  conducting  research  towards  integrating 
computational  fluid  dynamics  (CFD)  codes  with  computational  structural  mechanics 
(CSM)  codes.  This  will  give  the  Directorate  the  ability  to  develop  solutions  to  aeroelastic 
design  and  analysis  problems  using  state-of-the-art  computational  methods.  The 
integration  of  the  CFD  and  CSM  equations  is  done  cyclically.  The  fluid  dynamics 
equations  are  solved  time-accurately  for  a  single  time  step.  Next,  the  deflections  of  the 
body  are  computed  by  the  CSM  code  using  surface  pressures  computed  by  the  CFD  code. 
The  new  surface  deflections  are  then  used  to  update  the  aerodynamic  grid.  With  the 
updated  grid,  the  CFD  code  is  executed  for  another  time  step  and  the  process  repeats. 
This  paper  presents  an  algorithm  for  updating  the  grid  point  positions  of  the  aerodynamic 
grid.- 

An  algorithm  has  been  developed  for  deforming  aerodynamic  grids  that  is  based  on 
a  spring  analogy  (References  1  and  2).  This  method,  however,  is  seen  to  have  two 
limitations.  The  first  is  its  reliance  on  grid  connectivity.  The  computer  code  for  a 
structured  grid  would  be  different  from  that  for  an  unstructured  grid  and  the  code  for  a 
two-dimensional  grid  would  be  different  from  that  for  a  three-dimensional  grid.  The 
second  limitation  is  that  the  grid  point  deflections  are  a  function  of  local  grid  point 
spacing,  rather  than  body  geometry  and  global  distances.  Consequently,  the  method  best 
maintains  orthogonality  in  regions  of  high  grid  point  densities,  regardless  of  the  proximity 
to  the  surface. 

A  second  algorithm  has  been  developed  by  this  author  that  eliminates  the 
restrictions  of  the  spring  analogy  method  (Reference  3).  The  computer  code  can  be 
applied  to  a  wide  variety  of  grid  types  and  body  geometries  without  modification.  The 
algorithm  can  accurately  solve  for  large  deflections.  Also,  the  orthogonality  of  grid  lines 
at  a  surface  is  made  a  function  of  the  global  geometry:  the  deflection  at  a  given  point  in 
the  aerodynamic  grid  is  a  function  of  the  deflection  of  the  entire  aerodynamic  surface. 
However,  this  global  geometry  dependence  is  computationally  expensive.  For  three- 
dimensional  grids,  the  computational  cost  is  prohibitive.  Another  algorithm  is  needed 
which  is  both  computationally  efficient  and  able  to  produce  a  well-deformed  grid. 
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Presented  herein  is  a  third  grid  deformation  algorithm  that  is  based  on  a  single 
first-order  differential  equation.  The  smoothing  at  a  point  is  dependent  on  the  distance 
from  the  body.  The  equation  is  non-iterative  and  is  solved  by  marching  away  from  the 
aerodynamic  body.  The  method  is  computationally  inexpensive,  yet  maintains  the  quality 
of  the  original  undeformed  grid.  The  algorithm  is  dependent  on  the  grid  topology, 
however.  As  such,  the  grid  deformation  code  must  be  specialized  to  a  particular  grid 
topology. 


Grid  Deformation  Algorithm 


The  deformation  algorithm  solves  a  non-linear  relationship  between  known 
(independent)  boundary  deflections  and  unknown  (dependent)  interior  grid  point 
deflections.  The  deformation  algorithm  is  a  non-iterative  method  based  on  an  ad-hoc  first- 
order  differential  equation.  The  algorithm  is  applied  at  each  iteration  of  the  CFD-CSM 
loop. 


The  algorithm  can  be  divided  into  five  steps: 

1 .  If  desired,  transform  the  known  grid  point  deflections  to  deflections  with  respect  to  a 
body-fixed  coordinate  system. 

2.  Determine  appropriate  marching  directions. 

3.  Solve  for  the  unknown  grid  point  deflections. 

4.  If  necessary,  transform  the  grid  point  deflections  back  to  the  global  coordinate  system. 

5.  Update  the  grid  point  positions. 

Each  step  is  described  in  the  following  sub-sections. 

Preliminaries 

The  aerodynamic  surfaces  and  the  surrounding  grid  are  defined  in  terms  of  physical 
coordinates  x  .  To  facilitate  the  solution  of  the  CFD  equations,  the  grid  is  mapped  onto  a 
uniform  rectilinear  coordinate  system  ^ .  The  engineer  constructs  this  mapping  by  using 
an  analytic  or  computational  grid  generator.  Figure  1  shows  a  sample  mapping.  A  C-H 
grid  about  a  flat  plate  is  mapped  from  physical  coordinates  x  to  computational 
coordinates  ^  .  The  flat  plate,  shaded  gray,  maps  onto  a  portion  of  the  E,  face  of  the 
computational  domain.  Diagonal  hash  marks  identify  the  region  behind  the  flat  plate. 

A  grid  deformation  algorithm  must  maintain  the  quality  of  the  original  grid.  The 
following  conditions  should  be  satisfied  for  grid  quality  to  be  maintained: 

a)  No  negative  cell  volumes. 

b)  Grid  point  clustering  is  maintained. 

c)  Orthogonality  of  grid  lines  at  a  physical  surface  is  maintained. 

d)  Minimal  change  in  cell  volumes  near  a  physical  surface. 
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Figure  1.  Physical  to  Computational  Coordinate  Transformation 


Euler  Coordinate  Transformation 

The  first  step  of  the  grid  deformation  algorithm  is  to  transform  the  known  grid 
point  deflections  to  deflections  with  respect  to  a  body  fixed  coordinate  system.  This  is  an 
optional  step  that  helps  to  maintain  the  orthogonality  of  grid  lines  at  the  physical  surface, 
which  is  an  important  property  for  many  CFD  codes.  If  the  transformation  is  used,  the 
deflections  on  the  outer  boundary  will  not  be  zero  upon  transformation  back  to  the  global 
coordinate  system.  A  moving  outer  boundary  may  be  incompatible  with  physical 
constraints  or  may  require  special  boundary  conditions.  See  page  11  for  an  example  of  a 
grid  deformed  with  and  without  a  coordinate  transformation. 


An  Euler  transformation  is  used  to  transform  the  grid  point  deflections  to  the 
body-fixed  coordinate  system.  The  transformation  is  defined  as  three  consecutive 
rotations  with  respect  to  the  body-fixed  coordinate  system  about  a  center  point  defined 
in  the  global  coordinate  system.  The  rotation  angles  about  the  x, ,  .r, ,  and  Xj  body  axes 
are  6,,  0^,  and  63,  respectively,  and  have  corresponding  transformation  matrices  [©,], 
[Gt],  and  [03].  It  should  be  stressed  that  the  order  of  rotation  is  important.  The 
equation  for  the  transformed  deflections  of  points  on  the  physical  surface  is 


}  =  [0, I©: I©, ]{x"  +  u”  -  Xo } - {x;  - Xo } . 
The  product  of  the  three  transformation  matrices  is 


(1) 


[©.][©:  I©3]  = 


ce.ce, 

C0,s03  -sGiSe.ce, 

50,.S03  -t-C0,503C03 


-00,503 

C0,C03  +50,50,503 
50,003-00,50,503 


-50, 

-50,00, 

00,00, 


(2) 
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where  C  and  S  are  shorthand  for  the  cosine  and  sine  functions,  respectively.  If  the 
coordinate  transformation  is  not  applied,  then  Equation  1  simplifies  to 

{<}={<}  (3) 

since  the  rotation  angles  are  all  zero. 

Selection  of  Marching  Directions 

The  differential  equation  governing  the  grid  point  deflections,  described  in  the  next 
sub-section,  requires  an  initial  value.  Deflections  are  computed  along  grid  lines  by 
marching  away  from  the  initial  value.  However,  in  general,  there  is  not  a  single 
computational  coordinate  that  can  serve  as  the  marching  direction  for  the  entire 
aerodynamic  grid.  In  such  cases,  it  is  necessary  to  divide  the  aerodynamic  grid  into 
domains  within  each  of  which  there  is  a  consistent  marching  direction.  Often,  the  solution 
of  one  domain  will  be  used  as  an  initial  condition  for  a  second  domain. 

As  an  example,  consider  a  two-dimensional  C-grid  about  an  airfoil,  as  shown  in 
Figure  2.  The  chordwise  direction  is  and  the  normal  direction  is  Define  the  wake- 
cut  as  the  set  of  points  on  the  grid  line  that  connect  the  trailing  edge  with  the  outer 
boundary.  Deflections  on  the  airfoil  surface  are  known  from  the  CSM  code.  Deflections 
are  to  be  determined  for  the  points  in  the  wake-cut,  the  grid  interior,  and  the  outer 
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Figure  2.  C-Grid  About  an  Airfoil 
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boundaries.  The  wake-cut  points  are  chosen  to  form  one  domain  and  the  interior  points 
and  outer  boundary  points  a  second  domain.  The  deflections  in  the  wake-cut  domain  are 
obtained  by  using  the  deflection  of  the  trailing  edge  as  an  initial  condition  and  as  the 
marching  direction.  The  second  domain  uses  both  the  airfoil  and  wake-cut  grid  point 
deflections  as  initial  conditions  and  the  marching  direction  is  • 


Solving  for  Interior  Grid  Point  Deflections 


Let  4  2  be  the  marching  direction.  The  grid  point  deflections  in  terms  of  the  body- 
fixed  coordinate  system  are  determined  according  to  the  following  ad-hoc  non-linear  first- 
order  differential  equation: 


du' 


3^. 


-v(Jc,^) 


dx 

U 

(4) 


Each  term  is  described  below.  The  method  is  considered  ad-hoc  since  the  grid 
deformation  problem  is  not  based  on  any  natural  phenomenon  and  is  constructed  purely 
based  on  considerations  of  computational  accuracy  and  efficiency.  Equation  4  is  solved  by 
marching  along  grid  lines  of  constant  and  ^3 ,  with  varying  from  to  ^2^  •  Th® 
deflection  of  the  grid  point  serves  as  a  boundary  or  initial  condition.  The  ^2 min  point 
will  be  either  a  point  on  the  physical  surface  or  a  point  in  another  domain  for  which 
deflections  have  already  been  computed. 


Equation  4  produces  the  desired  behavior  in  the  region  away  from  the  physical 
surface:  the  deflection  u"'  and  the  gradient  of  the  deflection  BM'”/d^2  ^0  of  opposite 
sign  since  every  term  in  Equation  4  other  than  is  always  positive.  Also,  as  the 
deflection  u"'  goes  to  zero,  so  does  the  gradient  dM^/3^2  ■ 


Equation  4  can  be  written  as 


du' 


where  K  is  a  function  of  x  ,  ^  ,  and  U  : 

K(T,it/) 


-  =  -KM 

v(x,|) 

dx 

U 

9^2 

(5) 


(6) 
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If  K  is  replaced  by  a  constant.  Equation  5  becomes  linear  and  has  an  exponentially 
decaying  solution.  Thus,  it  is  clear  that  the  rate  of  decay  of  the  deflections  is  related  to  the 
function  k(x,^)  . 

The  spacing  function  v  is  chosen  to  give  the  desired  grid  behavior  near  the 
physical  surface:  surface  deflections  are  propagated  some  distance  into  the  interior  of  the 
aerodynamic  grid  with  little  smoothing.  The  function  v  is  plotted  in  Figure  3  and  is 
defined  as 


v(x,^)  = 


0 


/  \3 


V  'I  y 


f  / 

2- 


r.  +  -  r 


A 


V 


V  h  y 


r  <  To  -  r, 


ro  -  r,  <  r  <  To 


r^<r<r,+r, 


/■o  +  r,  <  r 


(7) 


Figure  3.  Grid  Spacing  Function 
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where  and  are  two  constants  chosen  by  the  user.  The  constant  is  a  measure  of 
the  distance  from  the  body  at  which  smoothing  becomes  significant.  The  constant  r,  is  a 
measure  of  how  rapidly  the  smoothing  becomes  significant.  No  smoothing  occurs  at 
distances  less  than  -  r,  away  from  the  surface.  The  function  r  is  defined  as  the 
distance  between  the  point  at  the  base  of  the  grid  line  and  the  interior  grid  point, 


and  is  as  shown  in  Figure  4. 

The  term  U  is  taken  to  be  slightly  greater  than  the  magnitude  of  the  maximum 
deflection  on  the  physical  surface.  The  magnitude  is  evaluated  according  to  the 

equation 


dx 

The  |3x/3^2|  U  terms  exist  in  Equation  4  to  help  prevent  the  creation  of  negative  ceil 
volumes  in  the  deformed  aerodynamic  grid. 


(9) 
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^2 


Figure  5.  Distances  Between  Consecutive  Grid  Points, 
Undeformed  and  Deformed 


Along  a  gridline,  two  consecutive  points  will  not  cross  if  the  following  condition  is 
satisfied: 


Consider  a  grid  line,  as  shown  in  Figure  5.  The  deflection  at  point  n  is  taken  to  be 
known  and  the  deflection  at  point  n  + 1  is  to  be  computed.  The  undeformed  distance 
between  the  two  points  is  \Ax\ .  The  deformed  distance  is  |Aic  +  Am|  .  The  magnitude  of 
Au  must  be  less  than  the  magnitude  of  Aic  for  the  g^d  line  to  remain  uncrossed. 
Letting  the  distance  between  the  two  points,  \Ax\ ,  go  to  zero,  the  condition  becomes  that 
of  Equation  10.  Equation  10  is  satisfied  by  Equation  4  since  Iw^/C/I  and  v  are  always 
less  than  one.  This  property  does  not  prevent  two  g^d  lines  from  crossing.  However, 
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grid  lines  should  not  cross  if  the  physical  surface  deforms  smoothly  and  if  deflections  do 
not  change  rapidly  over  the  surface. 

Inverse  Coordinate  Transformation 

With  all  of  the  grid  point  deflections  known  in  terms  of  the  body-fixed  coordinate 
system,  the  grid  point  deflections  in  terms  of  the  global  coordinate  system  can  be 
determined  by  applying  an  inverse  coordinate  transformation: 

(s; }  =  [0,  ]"  [©,  Ho ,  f  p:  + -  Jo }  -  {*:  -  ^0 }  ■  (11) 

Updating  Grid  Point  Positions 

With  the  grid  point  deflections  known  for  the  entire  grid,  the  grid  point  positions 
can  be  updated  according  to 

Deflections  for  the  (m+l)^^'  iteration  of  the  CFD-CSM  loop  are  referenced  to  the  grid  at 
the  iteration. 


Numerical  Implementation  of  the  Smoothing  Equation 


The  differential  equation  governing  the  grid  point  deflections.  Equation  4,  must  be 
solved  numerically.  The  partial  derivative  can  be  replaced  by  a  first-order  forward  finite 
difference.  Equation  4  then  becomes 


,  y  dx  Vm 


Since  the  grid  point  deflections  need  not  satisfy  any  physical  laws,  a  first-order  finite 
differencing  scheme  is  adequate. 
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Sample  Deformed  Grids 


This  section  gives  results  from  two  test  grids  deformed  with  the  method  presented 


in  this  paper.  The  first  is  a  two-dimensional 
0-grid  about  a  cylinder.  The  second  is  a 
three-dimensional  C-H  grid  about  a  swept 
wing. 

An  O-Grid  About  a  Cylinder 

Figure  6  shows  an  0-grid  about  a 
cylinder  of  unit  radius.  The  outer  boundary  is 
a  distance  of  5  radii  away  from  the  surface  of 
the  cylinder.  Grid  points  are  clustered  near 
the  surface  of  the  cylinder.  In  all  of  the 
deformed  grids  about  the  cylinder,  values  of 
/-Q  =  25  and  =  3.0  were  used  in  Equation 
4. 


Figure  6.  O-Grid  About  a  Cylinder 


Figure?.  Deformed  O-Grid  About  a  Figure  8.  Deformed  O-Grid  About  a 

Rotated  Cylinder  without  Coordinate  Rotated  Cylinder  with  Coordinate 

Transformation  Transformation 
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A  deformed  grid  resulting  from  a  30°  rotation  of  the  cylinder  is  shown  in  Figure  7. 
This  grid  was  obtained  without  applying  a  transformation  to  a  coordinate  system  fixed  to 
the  cylinder.  The  clustering  of  grid  points  near  the  surface  of  the  cylinder  has  been 
maintained.  However,  the  grid  lines  are  no  longer  normal  to  the  surface  of  the  cylinder. 


Figure  5  shows  the  effect  of  the  coordinate 
transformation  on  the  deformed  grid.  The 
entire  grid  is  simply  rotated  by  the  same 
amount  as  the  cylinder. 

Figure  9  shows  the  result  of  a 
translation  of  the  cylinder.  The  cylinder  is 
translated  half  a  radius  in  the  negative  x- 
direction  and  three-halves  a  radius  in  the  y- 
direction.  In  this  case,  both  clustering  and 
orthogonality  have  been  maintained.  Since 
the  body-fixed  coordinate  system  is  not 
rotated  with  respect  to  the  global  coordinate 
system,  the  transformation  would  have  no 
effect.  The  smooth  transition  between  the 
inner  and  outer  boundaries  of  the  grid  is 


Figure  9.  Deformed  O-Grid  About  a 
Translated  Cylinder  without  Coordinate 
Transformation 


Figure  10.  Deformed  O-Grid  About  a 
Translated  and  Rotated  Cylinder  without 
Coordinate  Transformation 
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Figure  11.  Deformed  O-Grid  About  a 
Translated  and  Rotated  Cylinder  with 
Coordinate  Transformation 


easily  visible  in  this  figure. 


Two  final  examples  of  deformed  grids  about  the  cylinder  are  shown  in  Figure  10 
and  Figure  11.  Here,  the  cylinder  is  both  translated  and  rotated.  The  grid  in  Figure  10  is 
produced  without  applying  a  coordinate  transformation.  The  grid  in  Figure  11  was 
produced  using  a  coordinate  transformation.  In  Figure  10,  the  outer  boundary  has  not 
been  affected,  while  in  Figure  1 1  the  outer  boundary  has  been  rotated  by  the  same  amount 
as  the  cylinder.  Both  have  maintained  grid  clustering  near  the  cylinder  surface. 
Orthogonality  was  only  maintained  when  a  coordinate  transformation  was  applied. 

These  grids  show  that  the  differential  equation  of  Equation  4  maintains  grid  quality 
under  translational  motion  of  the  physical  surface,  but  that  the  quality  of  the  grid 
deteriorates  under  rotational  motion.  The  coordinate  transformation  of  Equation  1 
improves  the  quality  of  the  deformed  grid  in  cases  with  rotational  motion,  but  it  does  so 


T 
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by  rotating  the  entire  grid,  which  moves  the  outer  boundaries. 


C-H  Grid  About  a  Swept  Wing 

Figure  12  shows  a  portion  of  a  C-H  grid  about  a  swept  wing.  The  wing  was 
deformed  by  applying  enough  of  the  first  bending  mode  to  deform  the  wing  tip  by  5%  of 
the  root  chord.  Values  of  =  2.0  and  r^  =  Of  were  used  in  Equation  4  and  no 
coordinate  transformation  was  applied.  For  this  grid,  there  was  no  consistent  normal 
direction  in  computational  space.  Therefore,  the  deformation  problem  was  divided  into 
three  domains.  First,  the  deflections  of  the  grid  points  outboard  of  the  wing  tip  and  in  the 
plane  of  the  wing  are  computed  by  marching  away  from  the  wing  tip  in  the  spanwise,  or 
^2,  direction.  Second,  the  deflections  in  the  wake  cut  are  determined  by  marching  in  the 
chordwise,  or  direction.  At  this  stage,  all  of  the  deflections  in  the  ^3  plane  are 
known.  Finally,  the  remaining  grid  point  deflections  are  determined  by  marching  in  the  ^3 
direction. 


Figure  13.  Deflection  Contours  of  a  Deformed  Grid  About  a  Swept  Wing 
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Figure  13  shows  contours  of  the  grid  point  deflections  in  the  ^3^,^  plane  of  the 
grid.  The  wing  is  outlined  with  a  heavy  border.  The  deflections  remain  approximately 
constant  for  distances  less  than  away  from  the  wing.  The  grid  deflections  then  smooth 
out  to  zero  over  a  distance  of  a  few  root  chords.  The  roughly  rectangular  shape  of  the 
contours  is  a  result  of  using  multiple  marching  directions. 
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Conclusions  and  Recommendations 


This  memorandum  presents  an  algorithm  for  deforming  an  aerodynamic  grid  given 
deflections  on  a  physical  boundary.  The  algorithm  has  been  implemented  in  FORTRAN 
and  has  been  successfully  tested  on  both  two-  and  three-dimensional  grids. 

The  method  is  capable  of  producing  reasonable  grid  deflections  in  a  reasonable 
amount  of  computer  time  for  a  range  of  practical  problems.  The  method  is  faster  than 
both  the  spring  analogy  and  the  global-geometry  methods.  It  is  applicable  to  a  variety  of 
physical  configurations,  including  a  wing  and  a  wing-body.  Also,  the  method  can  be 
applied  to  a  variety  of  grid  schemes,  such  as  C-H  grids  and  H-H  grids.  It  is  necessary, 
however,  to  code  each  configuration  separately. 

The  method  is  limited  to  reasonably  smooth  deflections  since  smoothing  is  only 
applied  in  the  direction  normal  to  the  surface.  A  discontinuity  in  the  surface  deflections, 
such  as  occurs  at  a  deflected  flap,  is  propagated  away  from  the  surface  without  smoothing. 
The  method  could  be  improved  by  adding  a  term  to  the  differential  equation  to  smooth 
such  discontinuities  in  surface  deflections. 

The  method,  as  presented  here,  cannot  enforce  grid  point  deflections  on  outer 
boundaries.  In  many  cases,  the  grid  point  deflections  will  damp  to  zero,  but  the  algorithm 
does  not  impose  such  a  constraint.  It  should  be  possible  to  enforce  deflections  on  multiple 
boundaries,  though,  by  applying  this  algorithm  iteratively,  alternately  marching  away  from 
the  physical  surface  and  the  outer  boundary.  This  would  make  the  method  applicable  to 
interior  flow  analyses  and  to  multi-element  airfoil  analyses. 
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